% Script to read 1 hour obs mag data, use 2000-2014 data
% compare it with world magnetic model and find rms

clear


S = dir('/Users/manojnair/data/obs_mag_data/BGS data/*.nc');

time_array_hour = (datenum(2000,1,1,0,30,0): (1/(24)): datenum(2014,12,31,23,30,0))';

f = figure(1);
set(f,'Position',[40          19        1557         937]);

% get the apex geomag lat
apex_obslat = load('/Users/manojnair/projects/wmm_error/obs_apex_lat_long.txt');

n_nodata_stn = 0;

for i = 1:length(S),
    
    %if abs(apex_obslat(i,4)) > 0 & abs(apex_obslat(i,4)) < 20 ,
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS data/' S(i).name];
    
    % get data from 2000 to 2014
    start_index = (datenum(1999,12,31) - datenum(1995,1,1) + 1) * 24;
    ndata = (datenum(2014,12,31) - datenum(2000,1,1) + 1) * 24;
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, start_index, ndata, 0);
    f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);
    h_data = sqrt(x_data.^2 + y_data.^2);
    i_data = atan(z_data./h_data) * (180/pi);
    if  strcmp(S(i).name(1:3),'CSY') ,
        d_data = unwrap(atan2(y_data,x_data)) * (180/pi);
        
    elseif strcmp(S(i).name(1:3),'DRV')
        d_data = unwrap(atan2(y_data,x_data)) * (180/pi) + 360;
        
    elseif strcmp(S(i).name(1:3),'DMC')
        d_data = unwrap(atan2(y_data,x_data)) * (180/pi);
    else,
        d_data = atan(y_data./x_data) * (180/pi);
        
    end;
    
    % store observatory metadata
    obs(i).lat = obj.geospatial_lat;
    obs(i).lon = obj.geospatial_lon;
    obs(i).name = S(i).name(1:3);
    
    % get the magnetic latitude
    path(path,'/Users/manojnair/projects/GIC/SHA_obs');
    [theta_gm, phi_gm] = gg2gm_new(90 - obs(i).lat,	obs(i).lon);
    obs(i).maglat = - (theta_gm - 90);
    
    
    L = isnan(x_data);
    
    if (sum(~L)) > 1000,
        %if 0,
        
        if  obj.geospatial_lon > 180,
            obj.geospatial_lon = obj.geospatial_lon - 360;
        end;
        
        icount = 1;
        % get world magnetic model data
        for decimal_year = 2000:0.1:2014.0,
            [wmm_xyz(icount,:) wmm_hc(icount), wmm_dc(icount), ...
                wmm_ic(icount), wmm_fc(icount)] = wrldmagm(0,obj.geospatial_lat,obj.geospatial_lon,decimal_year,...
                sprintf('%d',floor(decimal_year/5)*5));
            icount = icount + 1;
        end;
        
        
        
        % Process X data
        
        L = isnan(x_data);
        obs(i).x_ndata = sum(~L);
        data_array = x_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            x_spline=ppval(time_array_hour,sp);
            x_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            x_spline= sp(1) + sp(2) * time_array_hour;
            x_spline(L) = NaN;
        end;
        
        % get the stats
        wmm_x = interp1(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,1), time_array);
        obs(i).x_data_wmm_rms = rms(data_array - wmm_x);
        obs(i).x_spline_wmm_rms = rms(x_spline(~L) - wmm_x);
        obs(i).x_spline_resid_rms = rms(data_array-x_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_x);
        L = stats.w == 0;
        obs(i).x_data_wmm_rms_rob = rms(data_array(~L) - wmm_x(~L));
        
        % Process Y data
        
        
        L = isnan(y_data);
        obs(i).y_ndata = sum(~L);
        
        data_array = y_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            y_spline=ppval(time_array_hour,sp);
            y_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            y_spline= sp(1) + sp(2) * time_array_hour;
            y_spline(L) = NaN;
            
        end;
        
        wmm_y = interp1(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,2), time_array);
        obs(i).y_data_wmm_rms = rms(data_array - wmm_y);
        obs(i).y_spline_wmm_rms = rms(y_spline(~L) - wmm_y);
        obs(i).y_spline_resid_rms = rms(data_array-y_spline(~L));
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_y);
        L = stats.w == 0;
        obs(i).y_data_wmm_rms_rob = rms(data_array(~L) - wmm_y(~L));
        
        % Process Z data
        
        L = isnan(z_data);
        obs(i).z_ndata = sum(~L);
        
        data_array = z_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            z_spline=ppval(time_array_hour,sp);
            z_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            z_spline= sp(1) + sp(2) * time_array_hour;
            z_spline(L) = NaN;
            
        end;
        
        wmm_z = interp1(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,3), time_array);
        obs(i).z_data_wmm_rms = rms(data_array - wmm_z);
        obs(i).z_spline_wmm_rms = rms(z_spline(~L) - wmm_z);
        obs(i).z_spline_resid_rms = rms(data_array-z_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_z);
        L = stats.w == 0;
        obs(i).z_data_wmm_rms_rob = rms(data_array(~L) - wmm_z(~L));
        
        % Process F data
        
        L = isnan(f_data);
        data_array = f_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            f_spline=ppval(time_array_hour,sp);
            f_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            f_spline= sp(1) + sp(2) * time_array_hour;
            f_spline(L) = NaN;
        end;
        
        % get the stats
        wmm_f = interp1(datenum(2000:0.1:2014.0,0,0), wmm_fc, time_array);
        obs(i).f_data_wmm_rms = rms(data_array - wmm_f);
        obs(i).f_spline_wmm_rms = rms(f_spline(~L) - wmm_f);
        obs(i).f_spline_resid_rms = rms(data_array-f_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_f);
        L = stats.w == 0;
        obs(i).f_data_wmm_rms_rob = rms(data_array(~L) - wmm_f(~L));
        
        % Process DEC data
        
        L = isnan(d_data);
        data_array = d_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            d_spline=ppval(time_array_hour,sp);
            d_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            d_spline= sp(1) + sp(2) * time_array_hour;
            d_spline(L) = NaN;
        end;
        
        % get the stats
        wmm_d = interp1(datenum(2000:0.1:2014.0,0,0), wmm_dc, time_array);
        obs(i).d_data_wmm_rms = rms(data_array - wmm_d);
        obs(i).d_spline_wmm_rms = rms(d_spline(~L) - wmm_d);
        obs(i).d_spline_resid_rms = rms(data_array-d_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_d);
        L = stats.w == 0;
        obs(i).d_data_wmm_rms_rob = rms(data_array(~L) - wmm_d(~L));
        
        
        % Process INC data
        
        L = isnan(i_data);
        data_array = i_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            i_spline=ppval(time_array_hour,sp);
            i_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            i_spline= sp(1) + sp(2) * time_array_hour;
            i_spline(L) = NaN;
        end;
        
        % get the stats
        wmm_i = interp1(datenum(2000:0.1:2014.0,0,0), wmm_ic, time_array);
        obs(i).i_data_wmm_rms = rms(data_array - wmm_i);
        obs(i).i_spline_wmm_rms = rms(i_spline(~L) - wmm_i);
        obs(i).i_spline_resid_rms = rms(data_array-i_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_i);
        L = stats.w == 0;
        obs(i).i_data_wmm_rms_rob = rms(data_array(~L) - wmm_i(~L));
        
        
        % Process H data
        
        L = isnan(h_data);
        data_array = h_data(~L);
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array);
        if length(b1) > 1,
            sp=spline(b1,data_array'/spline(b1,eye(length(b1)),time_array'));
            h_spline=ppval(time_array_hour,sp);
            h_spline(L) = NaN;
            
        else,%data length <=1 year
            sp=robustfit(time_array,data_array);
            h_spline= sp(1) + sp(2) * time_array_hour;
            h_spline(L) = NaN;
        end;
        
        % get the stats
        wmm_h = interp1(datenum(2000:0.1:2014.0,0,0), wmm_hc, time_array);
        obs(i).h_data_wmm_rms = rms(data_array - wmm_h);
        obs(i).h_spline_wmm_rms = rms(h_spline(~L) - wmm_h);
        obs(i).h_spline_resid_rms = rms(data_array-h_spline(~L));
        
        [p,stats] = robustfit(1:length(data_array), data_array - wmm_h);
        L = stats.w == 0;
        obs(i).h_data_wmm_rms_rob = rms(data_array(~L) - wmm_h(~L));
        
        
        subplot(311);
        
        plot(time_array_hour,x_data);
        hold on;
        plot(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,1),'k*');
        plot(time_array_hour,x_spline,'r','LineWidth',2);
        legend('OBS','WMM', 'FIT');
        set(gca,'FontSize',16);
        
        ylabel('X nT');
        datetick('x');
        axis([datenum(2000,1,1), datenum(2015,1,1), -inf inf]);
        
        title(S(i).name(1:3));
        
        subplot(312);
        
        plot(time_array_hour,y_data);
        hold on;
        plot(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,2),'k*');
        plot(time_array_hour,y_spline,'r','LineWidth',2);
        
        set(gca,'FontSize',16);
        ylabel('Y nT');
        datetick('x');
        axis([datenum(2000,1,1), datenum(2015,1,1), -inf inf]);
        subplot(313);
        
        plot(time_array_hour,z_data);
        hold on;
        plot(datenum(2000:0.1:2014.0,0,0), wmm_xyz(:,3),'k*');
        plot(time_array_hour,z_spline,'r','LineWidth',2);
        
        set(gca,'FontSize',16);
        ylabel('Z nT');
        datetick('x');
        axis([datenum(2000,1,1), datenum(2015,1,1), -inf inf]);
        
        
        set(gcf, 'PaperPositionMode', 'auto');
        saveas(gcf,['/Users/manojnair/data/obs_mag_data/BGS data/plots/' S(i).name(1:3) 'WMM'],'epsc');
        clf
        
        
        %
        %                 subplot(311);
        %
        %                 plot(time_array_hour,f_data);
        %                 hold on;
        %                 plot(datenum(2000:0.1:2014.0,0,0), wmm_fc,'r.');
        %                 plot(time_array_hour,f_spline,'k');
        %                 legend('OBS','WMM', 'FIT');
        %                 set(gca,'FontSize',16);
        %
        %                 ylabel('F nT');
        %                 datetick('x');
        %                 axis([datenum(1995,1,1), datenum(2015,1,1), -inf inf]);
        %
        %                 title([S(i).name(1:3) sprintf('Apex Lat = %5.2f WMM-DATA RMS %5.2f DATA-SPLINE RMS %5.2f ', ...
        %                     apex_obslat(i,4),obs(i).f_data_wmm_rms, obs(i).f_spline_resid_rms)]);
        %
        %                 subplot(312);
        %
        %                 plot(time_array_hour,d_data);
        %                 hold on;
        %                 plot(datenum(2000:0.1:2014.0,0,0), wmm_dc,'r.');
        %                 plot(time_array_hour,d_spline,'k');
        %
        %                 set(gca,'FontSize',16);
        %                 ylabel('DEC degrees');
        %                 datetick('x');
        %                 axis([datenum(1995,1,1), datenum(2015,1,1), -inf inf]);
        %                  title([S(i).name(1:3) sprintf('Apex Lat = %5.2f WMM-DATA RMS %5.2f DATA-SPLINE RMS %5.2f ', ...
        %                     apex_obslat(i,4),obs(i).d_data_wmm_rms, obs(i).d_spline_resid_rms)]);
        %                 subplot(313);
        %
        %                 plot(time_array_hour,i_data);
        %                 hold on;
        %                 plot(datenum(2000:0.1:2014.0,0,0), wmm_ic,'r.');
        %                 plot(time_array_hour,i_spline,'k');
        %
        %                 set(gca,'FontSize',16);
        %                 ylabel('DIP degrees');
        %                 datetick('x');
        %                 axis([datenum(1995,1,1), datenum(2015,1,1), -inf inf]);
        %                  title([S(i).name(1:3) sprintf('Apex Lat = %5.2f WMM-DATA RMS %5.2f DATA-SPLINE RMS %5.2f ', ...
        %                     apex_obslat(i,4),obs(i).i_data_wmm_rms, obs(i).i_spline_resid_rms)]);
        %
        %
        %                 set(gcf, 'PaperPositionMode', 'auto');
        %                 saveas(gcf,['/Users/manojnair/projects/wmm_error/plots/' S(i).name(1:3) '_FDI'],'jpeg');
        
        
        
        
        display(S(i).name);
        %pause;
        %clf
        
        
    else,
        display('No data');
        display(S(i).name);
        n_nodata_stn = n_nodata_stn + 1;
        nodata(n_nodata_stn) = i;
    end;
    % end;
end;

save '/Users/manojnair/projects/wmm_error/observatoty_processed_data_Aug16' obs apex_obslat nodata

%% plot rms of crustal field


load  '/Users/manojnair/projects/wmm_error/observatoty_processed_data' obs apex_obslat nodata
maglat = abs(apex_obslat(:,4));
maglat(nodata) = [];

data_array_x = [obs.x_spline_wmm_rms];
data_array_y = [obs.y_spline_wmm_rms];
data_array_z = [obs.z_spline_wmm_rms];
data_array_f = [obs.f_spline_wmm_rms];
data_array_d = [obs.d_spline_wmm_rms].*wmmh;
data_array_i = [obs.i_spline_wmm_rms].*wmmf;
data_array_h = [obs.h_spline_wmm_rms];
nbin = 0;

for i = 5:10:85,
    nbin = nbin + 1;
    %     mid_lat = 45;
    %     lat_bin_rad = 45;
    %
    mid_lat  = i;
    lat_bin_rad = 5;
    % The exact calculation resulted in an RMS of 137.96 for 5 million points.
    
    L = maglat >= mid_lat-lat_bin_rad & maglat <= mid_lat+lat_bin_rad;
    lat_in_bin = maglat(L);
    

    
    data_in_bin_x = data_array_x(L);
    data_in_bin_y = data_array_y(L);
    data_in_bin_z = data_array_z(L);
    data_in_bin_f = data_array_f(L);
    data_in_bin_h = data_array_h(L);
    data_in_bin_d = data_array_d(L);
    data_in_bin_i = data_array_i(L);
    
    
    data_in_bin = sqrt((data_in_bin_x.^2 + data_in_bin_y.^2 + data_in_bin_z.^2));
    [p,stats] = robustfit(1:sum(L), data_in_bin);
    LL = stats.w < .90;
     %LL = zeros(size(data_in_bin_f));
   
    f_rms(nbin) = rms(data_in_bin_f(~LL));
    h_rms(nbin) = rms(data_in_bin_h(~LL));
    x_rms(nbin) = rms(data_in_bin_x(~LL));
    y_rms(nbin) = rms(data_in_bin_y(~LL));
    z_rms(nbin) = rms(data_in_bin_z(~LL));
    d_rms(nbin) = rms(data_in_bin_d(~LL));
    i_rms(nbin) = rms(data_in_bin_i(~LL));
    
    
    
    subplot(424);
    plot((lat_in_bin(~LL)), data_in_bin_x(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('X (nT)');
    hold on;
    subplot(425);
    plot((lat_in_bin(~LL)), data_in_bin_y(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('Y (nT)');
    hold on;

    subplot(426);
    plot((lat_in_bin(~LL)), data_in_bin_z(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('Z (nT)');
    hold on;
    subplot(423);
    plot((lat_in_bin(~LL)), data_in_bin_f(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('F (nT)');
    hold on;
    subplot(427);
    plot((lat_in_bin(~LL)), data_in_bin_h(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('H (nT)');
    hold on;
    subplot(422);
    plot((lat_in_bin(~LL)), data_in_bin_i(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('Inclination*F (Deg-nT)');
    hold on;
    subplot(421);
    plot((lat_in_bin(~LL)), data_in_bin_d(~LL),'r*');
    set(gca,'FontSize',16);
    xlabel('Corrected geomagnetic latitude (deg)')
    ylabel('Declination*H (Deg-nT)');
    hold on;
    clear L LL
end;


 subplot(424);
    plot(5:10:90, x_rms,'b.-','MarkerSize',20);
    subplot(425);
    plot(5:10:90, y_rms,'b.-','MarkerSize',20);

    subplot(426);
    plot(5:10:90, z_rms,'b.-','MarkerSize',20);
    subplot(423);
    plot(5:10:90, f_rms,'b.-','MarkerSize',20);
    subplot(427);
    plot(5:10:90, h_rms,'b.-','MarkerSize',20);
    subplot(422);
    plot(5:10:90, i_rms,'b.-','MarkerSize',20);
    subplot(421);
    plot(5:10:90, d_rms,'b.-','MarkerSize',20);


%% plot rms of the spline removed data
subplot(426);
plot(abs([maglat]), [obs.d_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.d_spline_resid_rms], 8);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));
set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('Declination (deg)');


% plot rms of the spline removed data
subplot(427);

plot(abs([maglat]), [obs.i_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.i_spline_resid_rms], 6);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('Inclination (deg)');



%
subplot(425);

plot(abs([maglat]), [obs.f_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.f_spline_resid_rms], 8);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('F  (nT)');

% plot rms of the spline removed data
subplot(421);

plot(abs([maglat]), [obs.x_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.x_spline_resid_rms], 6);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('X  (nT)');


%
subplot(422);

plot(abs([maglat]), [obs.y_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.y_spline_resid_rms], 8);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('Y  (nT)');

%
subplot(423);

plot(abs([maglat]), [obs.z_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.z_spline_resid_rms], 8);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('Z  (nT)');

%
subplot(424);

plot(abs([maglat]), [obs.h_spline_resid_rms],'r*');
hold on
p = polyfit(abs([maglat])', [obs.h_spline_resid_rms], 6);
plot(sort(abs([maglat])), polyval(p,sort(abs([maglat]))));

set(gca,'FontSize',16);
%title(sprintf('Global average Error is %5.2f%c from %d observatories',mean(robmean), char(176), ndata));
set(gca,'FontSize',16);
xlabel('Corrected geomagnetic latitude (deg)')
ylabel('H  (nT)');


%% Print out the crustal field error

load  '/Users/manojnair/projects/wmm_error/observatoty_processed_data' obs apex_obslat nodata
maglat = abs(apex_obslat(:,4));
maglat(nodata) = [];

data_array_x = [obs.x_spline_wmm_rms];
data_array_y = [obs.y_spline_wmm_rms];
data_array_z = [obs.z_spline_wmm_rms];
data_array_f = [obs.f_spline_wmm_rms];
data_array_dh = [obs.d_spline_wmm_rms].*wmmh;
data_array_d = [obs.d_spline_wmm_rms];
data_array_di = [obs.i_spline_wmm_rms].*wmmf;
data_array_i = [obs.i_spline_wmm_rms];

data_array_h = [obs.h_spline_wmm_rms];


fprintf('Mag Lat, 	 F  (nT),	  D(Deg),  D*BH (nT-Deg),	 X (nT),	 Y (nT),	 Z (nT),	 H (nT),	 I (Deg),	I*BF (nT-Deg), ndata/alldata\n');
for i = 5:10:85,
    
        mid_lat = 45;
        lat_bin_rad = 45;
    %
%      mid_lat  = i;
%     lat_bin_rad = 5;
    % The exact calculation resulted in an RMS of 137.96 for 5 million points.
    
    L = maglat >= mid_lat-lat_bin_rad & maglat <= mid_lat+lat_bin_rad;
    lat_in_bin = maglat(L);
    
    data_in_bin_x = data_array_x(L);
    data_in_bin_y = data_array_y(L);
    data_in_bin_z = data_array_z(L);
    data_in_bin_f = data_array_f(L);
    data_in_bin_h = data_array_h(L);
    data_in_bin_d = data_array_d(L);
    data_in_bin_dh = data_array_dh(L);
    data_in_bin_i = data_array_i(L);
    data_in_bin_di = data_array_di(L);
    
    % Robust rejection. To disable, commentout following three lines and
    % remove comment on LL = zeros(...
    data_in_bin = sqrt((data_in_bin_x.^2 + data_in_bin_y.^2 + data_in_bin_z.^2));
    [p,stats] = robustfit(1:sum(L), data_in_bin);
    LL = stats.w < .90;

    %LL = zeros(size(data_in_bin_f));

    
    f_rms = rms(data_in_bin_f(~LL));
    h_rms = rms(data_in_bin_h(~LL));
    x_rms = rms(data_in_bin_x(~LL));
    y_rms = rms(data_in_bin_y(~LL));
    z_rms = rms(data_in_bin_z(~LL));
    d_rms = rms(data_in_bin_d(~LL));
    dh_rms = rms(data_in_bin_dh(~LL));
     di_rmd = rms(data_in_bin_di(~LL));
    i_rms = rms(data_in_bin_i(~LL));
    
    fprintf('%d to %d, %5.0f,  %5.2f, %5.0f,%5.0f , %5.0f, %5.0f, %5.0f , %5.2f, %5.2f, %d/%d\n', ...
        mid_lat-lat_bin_rad, mid_lat+lat_bin_rad, f_rms, d_rms, dh_rms, x_rms, y_rms, z_rms, h_rms, i_rms, di_rmd, sum(~LL),length(LL));
    
    clear L LL
end;

% Grid Variation Error is defined as average declination error for
% geographic latitudes > |55|

obslat = abs([obs.lat]);
obslat(nodata) = [];
data_array_d = [obs.d_spline_wmm_rms];
L = obslat >= 55;
lat_in_bin = obslat(L);
    data_in_bin_x = data_array_x(L);
    data_in_bin_y = data_array_y(L);
    data_in_bin_z = data_array_z(L);
    data_in_bin_d = data_array_d(L);
% Robust rejection. To disable, commentout following three lines and
% remove comment on LL = zeros(...
data_in_bin = sqrt((data_in_bin_x.^2 + data_in_bin_y.^2 + data_in_bin_z.^2));
[p,stats] = robustfit(1:sum(L), data_in_bin);
LL = stats.w < .90;
dv_rms = rms(data_in_bin_d(~LL));

fprintf('GV RMS is %5.2f\n',dv_rms);

%% Print out the disturbnce field error
load  '/Users/manojnair/projects/wmm_error/observatoty_processed_data' obs apex_obslat nodata

maglat = abs(apex_obslat(:,4));
maglat(nodata) = [];

data_array_x = [obs.x_spline_resid_rms];
data_array_y = [obs.y_spline_resid_rms];
data_array_z = [obs.z_spline_resid_rms];
data_array_f = [obs.f_spline_resid_rms];
data_array_dh = [obs.d_spline_resid_rms].*wmmh;
data_array_d = [obs.d_spline_resid_rms];
data_array_i = [obs.i_spline_resid_rms].*wmmf;
data_array_h = [obs.h_spline_resid_rms];

fprintf('Mag Lat, 	 F  (nT),	D (Deg),  D*BH (nT-Deg),	 X (nT),	 Y (nT),	 Z (nT),	 H (nT),	 I (Deg),	 ndata\n');

for i = 5:10:85,
    
  mid_lat = 45;
  lat_bin_rad = 45;
    
%     mid_lat  = i;
%     lat_bin_rad = 5;
    
    L = maglat >= mid_lat-lat_bin_rad & maglat <= mid_lat+lat_bin_rad;
    lat_in_bin = maglat(L);
    
    data_in_bin_x = data_array_x(L);
    data_in_bin_y = data_array_y(L);
    data_in_bin_z = data_array_z(L);
    data_in_bin_f = data_array_f(L);
    data_in_bin_h = data_array_h(L);
    data_in_bin_d = data_array_d(L);
    data_in_bin_dh = data_array_dh(L);
    data_in_bin_i = data_array_i(L);
    
    LL = zeros(size(data_in_bin_f));
    
    f_rms = rms(data_in_bin_f(~LL));
    h_rms = rms(data_in_bin_h(~LL));
    x_rms = rms(data_in_bin_x(~LL));
    y_rms = rms(data_in_bin_y(~LL));
    z_rms = rms(data_in_bin_z(~LL));
    d_rms = rms(data_in_bin_d(~LL));
    dh_rms = rms(data_in_bin_dh(~LL));
    i_rms = rms(data_in_bin_i(~LL));
    
    fprintf('%d-%d, %5.0f, %5.2f,  %5.0f, %5.0f,%5.0f , %5.0f, %5.0f , %5.2f,  %d\n', ...
        mid_lat-lat_bin_rad,mid_lat+lat_bin_rad, f_rms, d_rms, dh_rms, x_rms, y_rms, z_rms, h_rms, i_rms, sum(~LL));
    
    clear L LL
end;


% Grid Variation Error is defined as average declination error for
% geographic latitudes > |55|

obslat = abs([obs.lat]);
obslat(nodata) = [];
L = obslat >= 55;
dv_rms = rms(data_array_d(L));

fprintf('GV RMS is %5.2f\n',dv_rms);


%% calculate H value at 2007

for i = 1:length(obs),
    
    obs_lon = obs(i).lon;
    if  obs_lon > 180,
        obs_lon = obs_lon - 360;
    end;
    [~, wmmh(i),~,~,wmmf(i)] = wrldmagm(0,obs(i).lat,obs_lon,2007.0,'2005');
    
end;


wmmh(nodata) = [];
wmmf(nodata) = [];


